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Abstract 

Construction of maximally localized Wannier functions (MLWFs) has been implemented within 
O the linear combination of pseudo-atomic orbital (LCPAO) method. Detailed analysis using ML- 

WFs is applied to three closely related materials, single benzene (Bz) molecule, organometallic 

E 

^ Vanadium-Bz infinite chain, and V2BZ3 sandwich cluster. Two important results come out from 

the present analysis: 1) for the infinite chain, the validity of the basic assumption in the mechanism 
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I of Kanamori and Terakura for the ferromagnetic (FM) state stability is confirmed; 2) for V2BZ3, an 
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I ^1 Bz is newly revealed: the on-site energy of p6 states of edge Bzs is higher than that of middle Bz, 

T— I which further reduces the FM stability of V2BZ3. 
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I. INTRODUCTION 



The electronic ground state of a periodic system is solved in a set of Bloch functions 
(BFs). They are eigenfucntions of both the Hamiltonian and lattice translation operators 
and characterized by two good quantum numbers n and k, the band index and crystal 
momentum, respectively. Though BFs are widely used in electronic structure calculations, 
they are difficult to be visualized due to their delocalized nature and hence do not offer an 
intuitive physical picture for chemical bonding and other local correlations. An alternative 
representation which can overcome these weaknesses is Wannier functions (WFs). Compared 
with BFs, WFs are localized in real space and constitute a description in terms of localized 
functions. The localization properties of WFs depend on the choice of phase factors of the 
BFs. For a group of isolated bands (isolated means this set of bands are connected among 
themselves by degeneracies in energy, but separated from others by finite energy gaps in 
the whole Brillouin zone), the degrees of freedom in phase factors of BFs are equivalent to 
unitary transformations among themselves at each k. Marzari and VanderbiltP developed 
a procedure which minimizes the spread of WFs (the second moment around their centers) 
by refining this degree of freedom. This procedure leads to WFs that are called maximally 
localized Wannier functions (MLWFs). If the bands of interest are not isolated and attached 
to, or cross with, other bands, a prescription for extracting the interested bands out of 
entangled bands is required. This disentanglement procedure was proposed by Souza, Mazari 
and VanderbiltP which servers as a pre-processing before refining the unitary transformations 
among the selected bands. Another method for constructing WFs with optimal localization 
properties is based on the A^th order muffin-tin-orbital (NMTO) method.!^ In this work, 
we only consider the former approach. 

MLWFs have stimulated intensive interests since it brings new hope to calculate several 
properties of materials which are quite hard to do within the representation of BFs. Since 
they are real in contrast to the complex BFs and well localized in real space, one can visu- 
alize them and gain intuitive physical insight into the nature of chemical bonding. It is also 
possible to extract some characteristic parameters such as the MLWFs' centers and spreads. 
The displacements of MLWFs' center is directly related with modern theory of polariza- 
tion. The hopping integrals among MLWFs from parameter-free first-principles calculations 
can be used to construct model Hamiltonians, or as a starting point for hDA+lP^^ or 
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LDA+DMFT^ calculations for strongly-correlated systems. Band structure interpolation 
based on Hamiltonian in WFs representation is quite efficient, which can be used for highly 
accurate integration in reciprocal spacd^such as that in calculating anomalous Hall effectP^*^ 
and electron-phonon couplingJ ^^ * ^^ ! MLWFs are also used for linear-scaling calculations for 
large systems. 

In this paper, we report briefly the implementation of MLWFs within OpenMX,'^! a 
ffist-principles electronic structure calculation software package, which is based on the lin- 
ear combination of pseudo-atomic orbital (LCPAO) basis functions and norm-conserving 
pseudopotentials within local density approximation (LDA) or generalized-gradient approx- 
imation (GGA). Since OpenMX is designed for large-scale ab initio calculations on parallel 
computers, our implementation is anticipated to allow a fast computation of MLWFs for a 
wide variety of materials such as biomaterials, carbon nanotubes, magnetic materials with 
different complex geometrical structures. 

In this work, one of the organometallic compounds vanadium-benzene sandwich-like 
complex,'^SIlSl VnBz.„+i, is reexamined by our newly generated MLWFs. This complex is 
one of the analogues of ferrocene, a prototype of metallocene. Experimentally, V„Bz„_|_i 
with n < 4 have been found to be one-dimensional cluster and have ferromagnetic (FM) 
ground state with the total magnetic moment increasing nearly linearly with the cluster 
size.'^ Inspired by these findings, several theoretical calculation j^^*^^' have been made on 
(VBz)„=oo, an ideal one-dimensional infinite chain. In those works, (VBz)„=oo is found to 
have highly stable FM ordering and shows half-metallic behavior. Double exchange was 
proposed to be the mechanism of FM ordering.^ On the other hand, by examining the elec- 
tronic structure from GGA+f/ calculations, we have proposed that the mechanism of FM 
stability should be that proposed by Kanamori and Terakura.l^ While compared with the 
infinite chain, finite V^Bz^+i clusters are found to have much weaker FM stability though 
the same mechanism is applicable. By using a simple tight-binding model, we have shown 
that absence of p-d hybridization in one side of edge Bz leads to magnetic polarization of 
edge Bz even for the AFM coupling of two V atoms, which reduces the total energy of AFM 
situation and thus reduces the FM stability energy largely.'^^'^ In the present work, a much 
more straightforward and quantitative analysis in the representations of MLWFs shows that 
there is another important role of the energy difference between the edge Bz and middle 
Bz, which further destabilizes the FM state against AFM one. By considering this, the 
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tight-binding model constructed now gives more consistent result with that of the direct 
first-principles calculation. 

In the following, we will describe technical issues regarding the construction of MLWF 
within the LCPAO method and then three closely related examples are studied, namely, Bz 
molecule, V-Bz infinite chain and V„Bz„+i {n=2) finite cluster, to demonstrate the successful 
applications of our implementation. Finally, we will conclude in Section IV. 

II. METHODOLOGY 

We will briefly introduce the theory of MLWFs. Only those aspects that are closely 
related with our LCPAO method will be described in detail. Details of general aspects can 
be found in the original papers Ref . [H and Ref . El Some other technical issues can be found 
in Ref. EU which introduces another implementation of constructing MLWF, wannier90.'^ 

A. Maximally localized Wannier functions 

The nth WF, n being the band index, localized in unit-cell at R is defined as Fourier 
transforms of an isolated band expressed by BF ^ as follows: 



where the integral is performed over the whole Brillouin zone (BZ). V is the volume of unit 
cell and N is the number of unit cells in the sample, e"*'^" '' is the undetermined phase 
factor, which brings indeterminacy of WF even transformed from a single isolated band. 
For a more general case with an isolated group of bands, e"*'^"-'' is generalized to a unitary 
transformation matrix U^^^: 



where is the total number of BFs in the isolated group of bands, the same as the number 
of WFs. In one single isolated band case, it has been proven that a suitable choice of the 
phase of e~^'^"''^ leads to WFs which are real and exponentially decaying in real space.'^ 
In multi-band case, the arbitrariness in the gauge transformation U^^^ can be exploited. 




(1) 
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According to Marzari and VanderbiltP among all of the arbitrary choices, a particular set 
will minimize the total spread of WFs, which is defined as 

n 

where and (r)„ are the expectation values of operators and r on the rath WF, respec- 
tively. Both expectations expressed in WFs can be transformed into those in BFs as shown 
by Blount.'^^ In practical calculation, a uniform k-grid is sampled to calculate the derivation 
of cell-periodic part of BFs in reciprocal space within finite difference approximation and the 
integral over k space is performed with summation over this grid. It is demonstrated that 
the dependence of Vt on the gauge transformation L/*-^^ is determined only by the so-called 
overlap integrals M^'^: 

M^^ = (V^„.,k|e-^^-'-|^n,k+b) = (n™,k(r)|M„,k+b(r)), (4) 

where Mm.klr) is the cell-periodic part of the Bloch states '?/'m,k = Wm,k(i")e~*'^''' and b is the 
vector connecting neighboring k-points in the regularly discretized mesh of k-points. 
is at the center of optimizing the spread of WFs since both the spread function itself and its 
gradient with respect to U^^'> are determined by it. Actual calculation of depends on 
the basis set used for electronic structure calculation and will be described in the following 
subsection. 



B. M^i in LCPAO method 

The wave-function within LCAPO method is defined as follow: 

1 ^ 

^..,k(r) = -7^Y. E - - R-P) (5) 

where 0ia(i" — Ti — Rp) is the pseudo-atomic orbital a centered on site in unit-cell Rp and 
^m\a is linear combination coefficients of them at k for band m. The overlap integral 
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matrix element is 



^mn = K,k(r)|ii„,k+b(r)) 
1 ^ 



(6) 



(0i,(r - - Rp)|e-^('-«-<'>''|0,-^(r - r,- - R,)) 



Defining that r' = r — — Rp, it becomes 



1 ^ 



(0,,,(r')|e-(r'+-'+^--^^)-^|(/),.;3(r' + r, - r,- + R, - R,)), 



(7) 



in which each term depends on only the relative position Rp — R^. Therefore, Eq. (7) can 
be written as 

N 

P ioi,jl3 (8) 

{4>ia{r')\e-''''-''\(pjp{r' + n- r,- + R^)). 

A uniform grid in real space is used to perform the integral. In practice, e"*'" '^ is expanded 
in terms of r' • b 

i=l i,j=i 

Since r' is an extended operator, we find that an expansion up to 4th order is needed to 
well conserve the unitary condition of M^^. A denser A;-space sampling will give smaller b, 
which helps to improve this conservation and lower order expansion can be used. 



C. Initial guess of MLWFs 

The minimization of the spread function begins with an initial guess for the target WFs. 
Following the approach proposed by Marzari and Vanderbilt, a set of trial functions \gn{r)), 
n e [1, Nu,], are taken as the initial guess of the N^, target MLWFs. In our LCPAO method, 
it is very convenient and natural to take the localized pseudo-atomic orbitals, which are the 
bases for expanding BFs, as the initial guesses. The center of each orbital can be put at any 
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place in the unit-cell and all the other characters for atomic orbital, such as the radial and 
angular functions, can be easily controlled by using suitable pseudo-atomic orbitals. The 
possible hybrids among the atomic orbitals are also available. 

D. General settings in OpenMX calculation 

Three closely related materials, single Bz molecule, V-Bz infinite chain and V2BZ3 cluster, 
are chosen to demonstrate our implementation. To calculate their electronic structures, 
PAOs are generated by a confinement potential scheme.^ For both hydrogen and carbon, the 
cutoff radius is 5.0 a.u. while it is 6.5 a.u. for vanadium. When generating pseudopotential, 
the semicore 3s and 3p states of V atom are included as valence states. The exchange 
correlation energy functional within GGJ^ is used for all the systems. Double- valence and 
polarization orbitals of each element are included as basis set: s2p2, s2p2dl and s2p2d2fl 
are used for H, C and V, respectively.l^Sl In the electronic structure calculation, the real-space 
grid techniqu^^ is used with an energy cutoff of 250 Ry in numerical integrations and in 
the solution of the Poisson equation. The GGA+U calculation is done with the approach 
proposed in Ref. EHl The geometrical structure of these materials are relaxed until the 
forces are less than 1.0 x 10""^ a.u. For molecular or cluster calculation, a supercell is used 
and the size is as large as 17 A to assure that the interaction between neighboring cells can 
be neglected. 

III. RESULTS AND DISCUSSIONS 
A. Benzene Molecule 

Benzene (Bz) molecule is firstly studied with only F-point sampled in BZ. Nine molecular 
orbitals (MO) around the HOMO and LUMO states are shown in Fig. 1 together with their 
eigenvalues and symmetries. Clearly HOMO-2, LUMO+2 and doubly degenerate HOMO, 
LUMO are composed of six pz orbitals on carbon atoms. All of the nine MOs are used to 
construct six MLWFs, i. e., Nwin=^ (Note that N^n denotes the number of band branches 
within the selected energy window as defined in Ref. |21) and Nyj=6. A physically intuitive 
initial guess for this set of target MLWFs are six orbitals on each carbon atom. The 
disentangling process proposed in Ref. |2] is used to select an optimized 6x6 subspace, 
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TABLE I: The first column shows the orbital energy of the pz like MLWF of a Bz molecule as an 
isolated one (first row). The rest of the columns are for the hopping integrals from wi to other 
five MLWFs. In the second row is for the V-Bz infinite chain. The numbers outside (inside) of the 
parenthesis are for spin up (spin down) channels. Energies are in eV. 





orbital energy 


hopping integrals between wi 


and Wi (i = 2 ~ 6) 




Wl 


W2,W6 W3,W5 


W4 


Bz 


-3.01 


-2.88 0.19 


-0.24 


(VBz)oo 


-4.23(-4.33) 


-2.64(-2.61) 0.08(0.08) 


-0.20(-0.19) 



which minimizes the gauge invariant part of the spread function, fi/. After that, a steepest- 
decent (SD) method is used to minimize the gauge-dependent part of the spread function to 
find the proper gauge transformation. With this initial guess, the total spread converged to 
IQ-^^ within 50 SD steps. The obtained six MLWFs are obviously identical to each other 
and have the spread of 0.943 A^, which is only slightly smaller than the initial spread of 
0.944 A2. MLWFs are real and have similar shape to the atomic pz orbital as shown in Fig. 
2. Although the initial guess of MLWFs are centered on each carbon atom, the converged 
MLWFs' centers are slightly shifted outward the gravity center of Bz by about 0.07 A along 
each nearby C-H bond direction. 

In Table. I, we listed the hopping integrals between these six MLWFs in the same unit 
cell. Wi means the MLWF centered around carbon atom i as shown in Fig. 2. The first 
column shows the orbital energy of the pz type MLWFs. The hopping integral from wi to its 
two nearest neighbors W2 and wq are -2.88 eV and those to further neighbors ws, and 
are 0.19 and -0.24 eV, respectively. The sign change in the hopping integral as the distance 
increases suggests that the WFs have oscillating tails to satisfy orthogonality relation. 

B. Ideal infinite (VBz) 

n=oa chain 

In our previous work,'^^*^ we have analyzed the mechanism of ferromagnetism stability 
in V-Bz infinite chain. Here we reexamine this system by using MLWFs. The spin polarized 
band structure FM V-Bz chain is shown in Fig. 3. Due to the symmetry, when V is 
sandwiched with Bz, its 3d orbitals will hybridize with Bz's HOMO-2, HOMO, LUMO to 
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LUMO+2 3,35eV(b2g) 



LUMO+1 3.30eV(a^g) 



LUMO -0.55 eV {e^J 



HOMO -5.S5 eV (e^^) 




HOMO-1 -7.81 eV (e^^) 




HOMO-2 eV (a^J 



FIG. 1: (color online). Molecular orbitals of Bz around the HOMO and LUMO states. The energy 
eigenvalue and symmetry of each orbital are shown. Doubly degenerate orbitals are shown together 
on the same level. 

form 3 types of bonds, namely a, vr and 6. Although HOMO-1 has also the 6 symmetry, its 
hybridization with dxy and d^i-y^ orbitals of V is negligible because the wave function of 
HOMO-1 is confined within Bz molecular plane. The 4s state of V is pushed to a high energy 
by about 5-6 eV with strong hybridization with HOMO-2. The original 4s electrons of V 
are transfered to the states near the Fermi level in Fig. 3. By including LUMO-l-2 orbital, 
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FIG. 2: (color online). The MLWF around the carbon atom 1 obtained for a Bz molecule is plotted 
with isovalue=0.1. All the six MLWFs are identical and their centers are shifted from nearby C 
atomic sites along each C-H bond by about 0.07 A. 

we construct eleven MLWFs (six Bz MOs and five V 3d-orbitals) from the eigenstates in the 
outer window from -8.4 eV to 6 eV as shown in Fig. 3. To demonstrate our implementation 
that can be applied to solid system, a 2 x 2 x 20 fc-grid is used. As in the Bz molecule 
case, the disentangling process is also necessary. We started from the initial guess of six 
orbitals on each carbon atom and five d orbitals on V atom. The tolerance for convergence 
is 10~^^ A^ for both gauge-invariant and -dependent parts of spread functions. Two schemes 
are used for optimizing the spread function. The SD method is used for first hundreds steps 
and it converg ed to IQ-^ A^. After that conjugate gradient (CG) method is adopted to 
continue the minimization, which takes about 40 steps to converge to 10~^^ A^. This hybrid 
scheme are found to be stabler than just using CG method in some cases, while faster than 
just using SD method. 

The quality of the obtained MLWFs can be seen from the interpolated band structure 
obtained by the band parameters for MLWFs shown in Fig. 3, which overlaps that of 
the original band at most k points. The obvious discrepancy happens in the a bands 
around the Fermi level originally composed of da and LUMO-l-1 states. This is because 
LUMO+1 state is discarded after disentangling. However, these two states have decreasing 
hybridization strength from F to X points, which leads to better reproduction of bands near 
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TABLE II: The spreads (in A^) of MLWFs of V-Bz infinite chain for GGA and GGA+U {U=3.0 
eV) calculations, pz means the Pz type MWLF, so as dz2, d^2_j^2, etc. $7/ is the gauge invariant 
part of the spread function. and 0,od represent the diagonal and off-diagonal contributions to 
the gauge dependent part of the spread function, respectively. 



spin 


Pz on Bz 


dz2 


dx^ — y^ :dxy 


dxZjdyz 








up(GGA) 


1.209 


1.155 


0.813 


1.052 


11.951 


0.001 


0.188 


down(GGA) 


1.186 


1.504 


0.870 


1.088 


12.264 


0.001 


0.273 


up(+C/) 


1.244 


1.080 


0.781 


1.096 


12.158 


0.001 


0.140 


down(-|-C/) 


1.179 


1.762 


0.935 


1.221 


12.814 


0.001 


0.341 



the X point .'^ The obtained MLWFs fulfill the requirement of real valuedness.'^ They are 
six Pz type orbitals centered nearby each carbon atom with a shift about 0.05 A along the 
C-H bond and five 3d-like orbitals centering on V atom site. The spreads for each MLWF 
and gauge decomposition are listed in Table II. The spreads of five d orbitals are classified 
into three categories corresponding to three types of bonds with Bz. The spread of dvr like 
MLWF is larger than that of d6 is due to the stronger hybridization with Bz's pz like MLWFs 
while the largest spread of da like MLWF is due to the da-dcr overlap along the chain. It 
is also noticed that the spreads of d orbitals in majority (up) spin channel are smaller than 
those in minority (down) one. As Bz has negative spin polarization against V, the spreads 
of spin down are also smaller than that of spin up, implying the occupied states may be 
more localized than unoccupied ones. 

The hopping integrals between these MLWFs can give more detailed information on 
bonding nature and physical mechanism in this system. The hopping integrals between 
six Pz orbitals are compared with those in an isolated Bz molecule in Table I. The on-site 
energy becomes -4.23 eV and -4.33 eV for spin up an spin down, respectively, which are 
lower than the one in an isolated Bz molecule due to the crystal field formed by ions. 
The hopping integrals between neighboring pz orbitals in one carbon ring are a little smaller 
than those in Bz molecule since the optimized C-C bond length in V-Bz chain is slightly 
larger than that in an isolated molecule The on-site energy of five d orbitals listed in Table 
III are also classified into three types corresponding to three types of bonds. The hopping 
integrals between d and Pz orbitals depend little on the spin polarization. Here only the 
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r X 



FIG. 3: (Color online) Band structure of FM V-Bz infinite chain. Lines are from ordinary first- 
principles calculations while symbols are those from Wannier interpolation. Red color is for spin 
up while green is for spin down. The component of interpolated bands at F point are labelled, too. 

hopping integrals between wi and d orbitals are listed. Those for other pz like MLWFs can 
be obtained by appropriate rotation of d orbitals. The pz-d hopping integrals decay quickly 
along the c-axis, the chain direction. For example, the Pz-d hopping integrals for the next 
nearest neighboring Bz are about 2 orders of magnitude smaller than those for the nearest 
neighbors and those for the third nearest neighbor can be ignored. This quick decaying 
property also indicates the good localization of obtained MLWFs. 

Similar calculation with f/=3.0 eV is performed to study how U influences the MLWFs. 
Compared with those in GGA calculation, +U makes the d-type (p^-type) MLWFs with 
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TABLE III: The onsite energies (second row) of d-type MLWFs and the Pz-d hopping integrals 
(third row) between one of the p^-type MLWF wi and d-type ones in FM V-Bz chain. In the 
fourth row, the Pz-d hopping integrals with p^ converted to molecular orbitals with a, vr and 6 
symmetries. The values outside (inside) the parentheses are spin up (spin down) channel and in 
unit of eV. The results are from the GGA calculations. 





d^2 










on-site 


-3.47(-1.99) 


-2.33(-1.35) 


-2.33(-1.35) 


-1.78(-1.21) 


-1.78(-1.21) 


Wl 


0.01(0.07) 


-0.77(-0.78) 


0.00(0.00) 


-0.96(-0.97) 


0.00(0.00) 


MO 


0.026(0.171) 


-1.34(-1.36) 


1.34(1.36) 


-1.66(-1.68) 


1.66(1.68) 



majority spin more localized (extended), while those with minority spin more extended 
(localized). The on-site energy of pz orbitals becomes -3.97 and -4.46 eV for spin up and 
spin down channels, respectively. The enhanced spin splitting of Pz orbitals is induced by the 
enhanced spin splitting in V d orbitals. The hopping integrals between different Pz orbitals 
keep nearly the same as those in GGA case. Similarly for V d orbitals, the spin splitting 
of the on-site energy is enhanced. For example for t/=3.0eV, the spin up (spin down) on- 
site energies of da, d5 and dvr states are -6.00 (-0.85), -3.91 (-0.37) and -1.38 (-0.27) eV, 
respectively. However, the pz-d hopping integrals have nearly no change. 

By using the MLWFs obtained so far, we reexamine quantitatively the tight-binding 
model of Ref. US] and Uni For this purpose, the atomic like MLWFs, Wi {i =1 to 6), 
are converted to MOs with symmetries of o", vr and 5. The essential point in Ref. [TS] 
are summarized below. HOMOs of Bz and dxz-, dj,^ orbitals having vr symmetry form quite 
strong bonds to hold the sandwich-like geometrical structure, while HOMO-2 and dz'^ orbitals 
having a symmetry couples weakly. The orbital with strong d^2 character is singly occupied 
and plays the role as a trigger of spin splitting. LUMOs and dx2_y2, d^y orbitals having 
6 symmetry form bonds to mediate magnetic coupling and they are mainly responsible for 
the stability of FM states. Our proposed mechanism for FM stability is described in Fig. 
4 following Kanamori and Terakura.l^ In this picture, the essential assumption is that the 
on-site energy of p5 states (LUMOs of Bz) should be between those of the spin majority 
and minority d states of 6 symmetry before switching on the p-d hybridization. Then, pS 
states could be negatively spin polarized due to the p-d hybridization leading to stability of 
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FIG. 4: (Color online) Schematic pictures describing the mechanism in which the p—d hybridization 
stabilizes the FM states. The DOS (a) without and (b) with p — d hybridization is plotted. The 
energy gain due to the charge transfer associated with the magnetic relaxation of p6 states is 
indicated. Before p — d hybridization, the on-site energies of p6 and d6 states are given in eV for 
both GGA and GGA+C/ calculations. See the text for more details. 

the FM configuration. The validity of the above assumption can be easily checked by the 
molecular orbitals formed with Wi {i =1 to 6). 

Similarly to a single Bz molecule case, diagonalizing 6x6 subspace Hamiltonian con- 
structed from the orbital energy and hopping integrals listed in Table I gives the eigen- 
energies and eigenstates of HOMO-2, HOMOs, LUMOs and LUMO+2. LUMOs, which 
play crucial roles in the FM stability have eigen-energies of -1.87 and -1.99 eV for spin up 
and spin down case, respectively. These two values are just sitting in-between the on-site 
energies of majority d6 state, -2.33 eV and minority one, -1.35 eV in GGA. This picture 
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is valid also in GGA+f/ calculation. In this case the spin up and spin down p6 states are 
located at -1.59 and -2.12 eV, respectively, which are well sitting between spin up d6 states 
at -3.91 eV and those of spin down at -0.37 eV. The p-d hopping integral between LUMOs 
and d6 {dx2-y2, dxy) orbitals can be easily obtained by using the hopping integrals between 
the atomic-like Pz MLWFs and d6 orbitals. For example, one of the LUMO states is given 
by 

\LUMO) = ^{\W2) - \ws) + \w,) - \we)). (10) 

Then in spin up case by using the values listed in Table III the hopping integral from this 
LUMO to dxy state is 

{dxy\H\LUMO) = 1.34. (11) 

Similarly, the hopping integrals between HOMOs and dvr orbitals and that between IIOMO-2 
and da orbital can also be calculated as listed in Table III. The bonding strength of different 
bonds can be seen from these hopping integrals. For cr and tt bonds, the values are around 
0.026(0.171) and 1.66(1.68) eV, respectively. The larger difference between spin up and spin 
down values in a bond case is a reflect of larger difference in spreads of d^2 like WF in 
spin up and down channels. Therefore, by the analysis using MLWFs, the physical pictures 
proposed in our previous work!^ are directly and quantitively confirmed. 

C. Finite V2BZ3 cluster 

We discussed the difference between infinite V-Bz chain and finite V„Bz„+i clusters in our 
previous two works.™Sl The reasons for nearly degenerate FM and AFM states in the finite 
V2BZ3 are extensively discussed. The importance of edge Bz is emphasized when compared 
with the infinite V-Bz chain. In the finite cluster, the edge Bzs having p-d hybridization 
only with one side of V atom can have magnetic relaxation even in the AFM configuration, 
which is not possible in the infinite chain. This makes the AFM state stabler and reduces 
the relative stability of FM state against AFM state. The p-d hybridization between V and 
edge Bz is stronger than that between V and middle Bz, which decreases the FM stability 
energy further. However, the present more quantitative analysis using MLWFs has revealed 
an additional source for making the FM and AFM states nearly degenerate in energy. 

The MLWFs for both FM and AFM V2BZ3 clusters have been constructed in the same 
way as above. We have found that MLWFs of edge Bz are different from those of middle 

15 



Bz in several ways. 1) The spreads of MLWFs from middle Bz are slightly larger than 
those from edge Bz by about 0.015 A^. 2) The hopping integrals from middle Bz to V d 
orbitals are slightly smaller than those from edge Bzs since V is closer to edge Bzs. This is 
consistent with our former analysis. 3) The on-site energy of LUMOs (p5 states) on edge Bz 
is about 0.78 and 0.88 eV higher than those on middle Bz for spin up and spin down states, 
respectively. The reason for the relatively higher p5 states on edge Bz is easily understood. 
Since the middle Bz is sandwiched by two V^^ ions while each of the edge Bzs have only 
one nearest V^"*" ion, the attractive electrostatic field on middle Bz is stronger than that 
on edge Bz and brings higher on-site energies for states on edge Bzs. While in our former 
work,'!^ we assumed that the on-site energies of p6 states are the same for both middle and 
edge Bzs. It is easy to check how the FM stability energy depends on the on-site energy 
difference between the edge Bzs and middle Bz by using the tight-binding model proposed in 
Ref. [15] In that model, only 6 bonds which are critical to the coupling of magnetic moments 
are considered. The parameters included in this model can be obtained directly from the 
construction of MLWFs. The on-site energy of p6 states on middle Bz, Sp, is found to be 
-2.34 eV (average value of those in spin up and spin down channels). Ed (on site energy of 
non spin polarized d6 state) and A (exchange parameter for d6 states) can be deduced from 
the orbital energy of spin polarized d6 orbitals. They are found to be -2.13 and 1.18 eV, 
respectively. Parameters for describing hybridizations between d6 orbitals and p5 states on 
edge Bzs and middle Bz are t' and t with the values of 1.44 and 1.28 eV, respectively. By 
using these parameters, it is shown in Fig. 5 that the FM stability energy decreases as the 
energy difference in p(5 states on edge Bzs and middle Bz increases. AFM state becomes 
stabler than FM state (FM stability energy becomes negative) when p6 states on edge Bz 
is about 0.67 eV higher than those on middle Bz. This critical value is a little bit different 
from the value, around 0.83 eV, obtained directly from MLWFs since in this tight-binding 
model other contributions, such as those from a bonds, are neglected. Therefore, we think 
that the difference between edge Bzs and middle Bz in p(5 states is another important source 
for the much reduced FM stability in finite clusters. 
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FIG. 5: (Color online) FM stability energy (total energy of AFM state minus that of FM state, in 
eV) depends on the energy difference in p6 states (Asp) on edge Bzs and middle Bz. The points 
where the FM stability energy equals 0.0 eV and the on-site energy difference being 0.83 eV are 
indicated by "+" symbols. 

IV. CONCLUSION 

We have implemented the construction of maximally localized Wannier funcitons in the 
formalism of linear combination of pseudo-atomic orbitals first-principles calculations. The 
implementation is demonstrated to be applicable to both solid and molecular systems. The 
disentangling procedure works well for metallic cases also. The analysis of hopping integrals 
obtained from MLWFs for V-Bz organometallic complex indicates that it is a proper way 
to find the tight-binding parameters from parameter-free ab initio calculations. MLWFs 
provide useful information for the study of the bonding nature and physical mechanism in 
materials. As demonstrated in V-Bz complexes, a new physical origin, the role of orbital 
energy difference between edge Bz and middle Bz neglected before, is naturally revealed by 
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analysis with MLWFs. 
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